syms u1 u2 u3 u4 u5 u6 mu
x1=u4;
x2=u5;
x3=u6;
x4=2*u5+u1-mu*(u1+mu-1)/(u2^2+u3^2+(u1+mu-1)^2)^(3/2)-(1-mu)*(u1+mu)/(u2^2+u3^2+(u1+mu)^2)^(3/2);
x5=-2*u4+u2-mu*u2/(u2^2+u3^2+(u1+mu-1)^2)^(3/2)-(1-mu)*u2/(u2^2+u3^2+(u1+mu)^2)^(3/2);
x6=-mu*u3/(u2^2+u3^2+(u1+mu-1)^2)^(3/2)-(1-mu)*u3/(u2^2+u3^2+(u1+mu)^2)^(3/2);
J=jacobian([x1;x2;x3;x4;x5;x6],[u1 u2 u3 u4 u5 u6]);
J